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Abstract: We present a new flexible, fast and accurate way to implement massive neutri- 
nos, warm dark matter and any other non-cold dark matter relics in Boltzmann codes. For 
whatever analytical or numerical form of the phase-space distribution function, the optimal 
sampling in momentum space compatible with a given level of accuracy is automatically 
found by comparing quadrature methods. The perturbation integration is made even faster 
by switching to an approximate viscous fluid description inside the Hubble radius, which 
differs from previous approximations discussed in the literature. When adding one massive 
neutrino to the minimal cosmological model, CLASS becomes just 1.5 times slower, instead 
of about 5 times in other codes (for fixed accuracy requirements) . We illustrate the flexibil- 
ity of our approach by considering a few examples of standard or non-standard neutrinos, 
as well as warm dark matter models. 
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1. Introduction 

The inclusion of massive, non-cold relics in a Boltzmann code is complicated by the fact 
that it is necessary to evolve the perturbation of the distribution function on a momentum 
grid. A grid size of N points together with L terms in the expansion of the perturbation 
leads to N ■ L added equations to the system. In public Boltzmann codes like CMBFAST Q, 
CAMB H and CMBEASY ||, distributions are sampled evenly with fixed step size and maxi- 
mum momentum, adapted to the case of a Fermi-Dirac shaped distribution function f(p). 
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Moreover, the analytic expression for f(p) is hard-coded in many places in those codes, 
and implicitly assumed e.g. in the mass to density relation, so that exploring other models 
like neutrinos with chemical potentials and flavour oscillations, neutrinos with non-thermal 
corrections, extra sterile neutrinos or any kind of warm dark matter candidate requires 
non-trivial changes to these codes. 

We present here the way in which generic Non-Cold Dark Matter (NCDM) relics are 
implemented in the new Boltzmann code CLASS 1 (Cosmic Linear Anisotropy Solving Sys- 
tem), already presented in a series of companion papers ||, g|. In order to ensure a 
complete flexibility, CLASS assumes an arbitrary number of NCDM species, each with an 
arbitrary distribution function fi(p). For each species, this function can be passed by the 
user under some (arbitrarily complicated) analytic form in a unique place in the code, or in 
a file in the case of non-trivial scenarios that requires a numerical simulation of the freeze- 
out process. All other steps (finding a mass-density relation, optimising the momentum 
sampling and computing the derivative of fi(p)) are done automatically in order to ensure 
maximum flexibility. 

In Sec. ^, we present an automatic quadrature method comparison scheme which allows 
CLASS to find an optimal momentum sampling, given fi(p) and some accuracy requirement, 
and in Sec. ||, we devise a new approximation scheme allowing us to drastically reduce the 
computational time for wavelengths inside the Hubble radius. Finally, in Sec. || and ||, 
we illustrate these methods with several examples based on standard and non-standard 
massive neutrinos, and different types of warm dark matter candidates. 



2. Optimal momentum sampling 



The formalism describing the evolution of any NCDM species is given by the massive 
neutrino equations of Ma & Bertschinger . We will follow the notations from this paper 
closely, with the exceptions 



Qm £mb (2,2 m2 i 2 

i=r ' e =r = \ q+a r2 — . 

^ncdm,0 ^ncdm,!} \ - t ncdm,0 / 



where T nc6m ^ is the temperature of the non-cold relic today, in the case of a thermal relic. 
If the relic is non-thermal, T ncdmi o is just a scale of the typical physical momentum of the 
particles today. Note that the perturbation equations Eq. (2.4) are still the same as in Q, 
since they depend only on the ratio q/e which is not affected by this rescaling. 



2.1 Perturbations on a grid 

We are not interested in the individual momentum components of the perturbation, 
but only in the perturbed energy density, pressure, energy flux and shear stress of each 

Available at http://class-code.net. This paper is based on version vl.l of the code. 
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NCDM species, which are integrals over ^ Q: 



^Pncdin = 47T 

(Pncdm Pncdia) ^ncdm 
(Pncdm Pncdm) "ncdm 



1 ncdm,0 

a 

4vr / T ncdm i0 



47r/c 



ncdra.O 



8vr / T T 



- ncdm,0 



f (q)dqq 2 e^ c 



fo(q)dq—ty , 



Mq)dqq 3 *i 



f (q)dq—y 2 - 



(2.2a) 
(2.2b) 
(2.2c) 
(2.2d) 



3 \ a 

In the rest of the article, we will omit all ncdm subscripts, and dots will denote derivatives 
with respect to conformal time, r. 

Note that *$>o an d ^i are gauge-dependent quantities, while higher momenta are not. 
The gauge transformation can be derived from the corresponding gauge transformation of 
the integrated quantities. The relation between *i in the conformal Newtonian gauge and 
in the synchronous one reads: 

lq 



*l,Con. = *l,Syn. + ok 



q 6 € 



(2.3) 



with a = (h + Qf])/(2k 2 ), where h and r\ are the usual scalar metric perturbations in the 
synchronous gauge. In the rest of this paper, we will work exclusively in the synchronous 
gauge. The evolution of the Vl^'s are governed by the Boltzmann equation as described 
in P], and leads to the following system of equations: 



qk „ h d In fn 
— _| ±1. 

e 6 ding 



qk 



#i = l_(tf -2tf 2 ), 
3e 



5e 



(2*i - 3* 



h 2rj\ d In / 



«>3 



V 15 5 I 



ding 



(2/ + 1) e 

We can write the homogeneous part of this set of equations as 



* = — = a (t) A^, 



(2.4a) 
(2.4b) 
(2.4c) 

(2.4d) 
(2.5) 



where A is given by 



A 



21+1 



l+l 
21+1 



(2.6) 
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The solution can be written in terms of the matrix exponential, 

^(T) = e f n dT ' a{T ' )A ^(n) (2.7) 
= Ue^ dT ' aiT ' )D U- l ^(r i ), (2.8) 

where A has been diagonalised such that A = UDU^ 1 and D is a diagonal matrix of 
eigenvalues of A. The largest eigenvalue of A (using the complex norm) goes toward ±i for 
l m£H — > oo, so the largest frequency oscillation in the system is 

r ( m 2 

Wmax^^y dr' [l + -p-a {t'YJ . (2.9) 
2.2 Quadrature strategy 

There is no coupling between the momentum bins, so our only concern is to perform the 
indefinite integrals numerically with sufficient accuracy while using the fewest possible 
points. We are interested in the integrals in Eq. fl2.2j), which are all on the form 

/•oo 

1= / dqf (q)g(q), (2.10) 
J 

where fo(q) is the phase space distribution and g(q) is some function of q. We will assume 
that g{q) is reasonably well described by a polynomial in q, which we checked explicitly for 
the functions in Eq. fl2.2|) . Under this assumption, we can determine the accuracy of any 
quadrature rule on I by performing the integral 

poo 

J= / dqf (q)t(q), (2.11) 
J 

where t{q) is a test function. Given a set of different quadrature rules for performing the 
integral X, the idea is to choose the rule which can compute J to the required accuracy 
tol_ncdm using the fewest possible points. 

We define a quadrature rule on I to be a set of weights W{ and a set of nodes qi, such 

that 

n 

Z~£Vi<7( fi ). (2.12) 

Note that the distribution function itself has been absorbed into the weights. The optimal 
quadrature rule will depend on both the distribution fo(q) and the accuracy requirement 
tol_ncdm, but the specific method used for obtaining the rule is decoupled from the rest 
of the code; the output is just two lists of n points, {qi} and {VFj}. CLASS tries up to 
three different methods for obtaining the most optimal quadrature rule, each with its own 
strength and weaknesses. These are Gauss-Laguerre quadrature, adaptive Gauss-Kronrod 
quadrature and a combined scheme. We will now discuss each of them. 
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2.3 Gauss-Laguerre quadrature 

Most of the time, the distribution function will be close to a Fermi-Dirac distribution, and 
the integrand is exponentially decaying with q. The Gauss-Laguerre quadrature formula 
is well suited for exponentially decaying integrands on the interval (0;oo), so this is an 
obvious choice. The rule is 

/ dqe- q h(q)^y j w l h(q l ), (2.13) 
J ° i=l 

where the nodes qi are the roots of L n , the Laguerre polynomial of degree n and the weights 
can be calculated from the formula 

Wi = _ (2.14) 

(n + l) 2 [L n+1 ( qi )} 2 

If we put h(q) = e q fo(q)g(q) we obtain the rule 

Wi = Wie«f (ft) . (2.15) 

This rule will be very effective when the ratio fo/e~ 9 is well described by a polynomial, 
but it will converge very slowly if this is not the case. 

2.4 Adaptive sampling 

When an integrand has structure on scales smaller than the integration interval, an adaptive 
integration scheme is often the best choice, since it will subdivide the interval until it 
resolves the structure and reach the required accuracy. We will use the 15 point Gauss- 
Kronrod quadrature formula as a basis for our adaptive integrator; 7 of the 15 points can 
be used to obtain a Gauss quadrature estimate of the integral, and the error estimate on 
the 15 point formula is then err est . = 200|G7 - K15\ 15 . 

The Gauss-Kronrod formula is defined on the open interval (—1,1), but it can be 
rescaled to work on an arbitrary open interval (a, b). We transform the indefinite integral 
into a definite integral by the substitution x = (q + l) -1 : 

J" dqf (q) = -£ dx^-f (q (x)) = £ dxx~ 2 f (q (x)) . (2.16) 

This integral can then be solved by the adaptive integrator. If the tolerance requirement is 
not met using the first 15 points, the interval is divided in two and the quadrature method 
is called recursively on each subinterval. 

This method is very efficient when the integrand is smooth. For practical purposes, 
this will be the case unless the phase-space distribution is read from a file with sparse 
sampling: in this case, the code must interpolate or extrapolate the file values in order to 
cover the whole momentum range, and the next method may be more efficient. 
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2.5 Integration over tabulated distributions 

If some distribution function is not known analytically, but only on a finitely sampled grid 
on (g m i n , (/max)) we have to interpolate the distribution function within the interval, and we 
have to extrapolate the behaviour outside the interval. Inside the interval we use a spline 
interpolation, while we assume f(q < q m i n ) = f(q m in) close to zero. For the tail, we assume 
the form f(q) = ae~^ q . Requiring the function and its first derivative to be continuous at 
the point q = g max leads to the following equations for a and (3: 



a = /(( / max)e /3<?max , (2.17) 

(2.18) 



P = -/(<7max) 1 ^ 



<2 — <?max 



In the combined scheme we use the 4 point Gauss-Legendre method on the interval (0, <7min)> 
adaptive Gauss-Kronrod quadrature on (q m in, Qmax.) and the 6 point Gauss-Laguerre rule 
on the tail (q ma , x , oo) 2 . This scheme works well when the integrand is interpolated from 
tabulated points. 

2.6 Implementation in CLASS 

When CLASS initialises the background structure, it will find optimal momentum samplings 
for each of the species. More specifically, we start by computing the integral of the distribu- 
tion function multiplied by the test function at high accuracy, which gives a reference value 
which can be used for comparison. It also creates a binary tree of refinements, from which 
we can extract integrals at various levels, where level 1 is the best estimate. We choose the 
highest possible level which results in an error which is less than the input tolerance, and 
we extract the nodes and weights from that level. 

The code will now search for the lowest number of nodes required for computing the 
integral with the desired accuracy using Gauss-Laguerre quadrature. The most efficient 
method, the method using the lowest number of points, is then chosen. For a distribu- 
tion not departing too much from a Fermi-Dirac one, this will usually be Gauss-Laguerre 
quadrature. 

The scheme suggested here has the benefit, that there is just one tolerance parameter 
directly related to how well the integral is approximated, independently of the distribution 
function. However, for this to be exactly true, we require the test function to be a suffi- 
ciently realistic representation of q n ^>i for n = 2, 3, 4 and I = 0, 1, 2 for the perturbations. 
We have checked this using different test functions, but in the end we found the polynomial 
f(q) = a,2q 2 + a^q 3 + a^ 4 to be adequate. The coefficients were chosen such that 

oo n 

dq— = 1. (2.19) 

o efi + 1 V ' 

When the phase-space distribution function is passed in the form of a file with tabulated 
(<7j, fj) values, the code compares the three previous methods (still with a common tolerance 



This version of the rule is obtained by a simple substitution. 
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parameter) and keeps the best one, which is usually the third one in the case of a poor 
sampling of the function, or one of the other two in the opposite case. 

Note that higher accuracy is needed for integrating background quantities (density, 
pressure, etc.) than perturbed quantities (the VlZ/'s). On the other hand, the code spends 
a negligible time in the computation of the former, while reducing the number of sampling 
points for perturbations is crucial for reducing the total computing time. Hence, CLASS calls 
the quadrature optimisation algorithm twice for each NCDM species, with two different 
accuracy parameters. The background tolerance is set to a smaller value leading to a finer 
sampling. 

We conclude this section by noticing that this whole process sounds very sophisticated, 
but requires a negligible computing time in CLASS. What really matters is to reduce the 
number of discrete momenta in the perturbation equations, and this is indeed accomplished 
thanks to the previous steps (as we shall see in Sec. |^). 



3. Sub-Hubble Approximation 
3.1 Fluid approximation 

Various kinds of approximations for massive neutrino perturbations have been discussed 



in the past |@, ||, 10]. The approximation discussed here is different and consists in an 
extension of the Ultra-relativistic Fluid Approximation presented in Q , applying only to the 
regime in which a given mode has entered the Hubble radius. The idea is that after Hubble 
crossing, there is an effective decoupling between high multipoles (for which power transfers 
from smaller Us to higher /'s, according to the free-streaming limit) and low multipoles 
(just sourced by metric perturbation). Hence, when kr exceeds some threshold, we can 
reduce the maximum number of multipoles from some high l max down to Z max = 2. We 
showed in [|| that this Ultra-relativistic Fluid Approximation (UFA) allows simultaneously 
to save computing time (by reducing the number of equations) and to increase precision 
(by avoiding artificial reflection of power at some large cut-off value £ m ax)- 

In the case of massive neutrinos, we expect the same arguments to hold in the relativis- 
tic regime, while in the non-relativistic limit all multipoles with I > 1 decay and the species 
behave more and more like a pressureless fluid. Hence, some kind of fluid approximation 
is expected to give good results in all cases. 

We write the continuity equation and the Euler equation in the usual way. In the 
synchronous gauge we have 

5 = - (1 + w) ( 9 + ~ ] - 3- (c| yn . - w) 5, (3.1a) 



= --{1-^6 + ^-^5-^0. (3.1b) 
a y 1 + w 

Here, c 2 g is the adiabatic sound speed, and c| yn = |jj is the effective sound speed squared 
in the synchronous gauge. The latter can be related to the physical sound speed defined 



7 



in the gauge comoving with the fluid, that we denote c e g. The above equations can then 
be written as: 



(i + w)\e + 



h 



W 



)« + 9(£) (l + w)(cl s -4) 



k 2 ' 



(3.2a) 



a 
a 



i-3c 2 eS )e + T 



-cff 



+ W 



■k 2 5 



k 2 a. 



(3.2b) 



Later on, we will close the system by an evolution equation for the shear a, but first we will 
discuss how to calculate the adiabatic sound speed and how to approximate the effective 
sound speed c 2 ^. 

3.2 Sound speeds 

The adiabatic sound speed can be expressed as 

-l 



V V P 
— = w— — 

P V\P 



w 



3(l + u;) 



p a 
-w- — 
V \a 



-i 



1 



sa + w) 



(3.3) 



where the quantity p (called the pseudo-pressure inside CLASS) is a higher moment pressure 
defined by 



47T 

p = — a 
F 3 



fo(q)dq—. 



(3.4) 



With this formulation, we can compute the adiabatic sound speed in a stable and accurate 
way, without needing to evaluate the time-derivative of the background pressure p. When 
the ncdm species is no longer relativistic, its pressure perturbation bp defined in Eq. ( |2.2b| ) is 
an independent quantity. Since we do not have an evolution equation for 5p, we approximate 
c eff by c g ■ This approximation is sometimes as much as a factor 2 wrong as shown on Fig. [I]. 




10 100 1000 10000 

conformal time (Mpc) 




10 100 1000 

conformal time (Mpc) 



Figure 1: Effective sound speed squared c^ ff for a mass of in = 2.0 eV. Left panel: The effective 
sound speed plotted together with the adiabatic sound speed squared c 2 g and the equation of state 



parameter w. 



In the relativistic and the non-relativistic limit we have 



c as expected, but 



the behaviour of c^ ff in between the two limits are non-trivial. Right panel: The ratios c\ s j cr g and 
c^ s /w. One can see that c 2 is a better approximation to c 2 s than w, but neither catches the full 



evolution. 
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3.3 Evolution equation for the shear 



Given an ansatz for ^3, we can derive a formally correct evolution equation for the shear. 
We follow Ma and Bertschinger, and close the system using their suggested recurrence 
relation for massive neutrinos 3 . The truncation law presented in Ma and Bertschinger is 
valid for Z max > 3: in this case, all quantities are gauge-invariant. When writing the same 
ansatz for Z max = 3, we have to face the issue of the gauge dependence of ^1. Asssuming 
that the truncation law holds for gauge-invariant quantities, one obtains in the synchronous 
gauge: 

5e / 
* 3 « ^—^2 - ( *i + Oik 



qkr 



e 1 q 
q 6 € 



(3.5) 



Throughout this subsection and the next one, one can recover Newtonian gauge equations 
by simply taking a = 0. We now differentiate equation fl2.2d| ): 



1 



57T 



p + p 3 



fo(q)dqq 4 



d_ 

dr 



^2 
€ 



(3.6) 



We can compute the right-hand side using Eq. ( |2.4c| ) and replace ^3 with its approximate 
expression from ( |3.5D . After carrying out integrals over momentum, one gets: 



a 



a 

+ - 

a 



3"^ 



1 E 
3 a 



2 

* + 3 



O + at 



3+ P - 
1 + w \ p 



where we have borrowed the notation 



f'°° a 2 
(p + p)Q = A7Tka' 4 f (q)dqq 3 ^i, 
Jo e 



(p + p)S = ^-a 



roo 4 2 

Jo € e 



(3.7) 

(3.8) 
(3.9) 



from flCf| . Erom the definition it is clear that — ?■ # and £ — > a in the relativistic limit, 
and that O and S become suppressed in the non-relativistic regime compared to 9 and 
o. Our differential equation for a differs from its Newtonian gauge counterpart in |]10|| , 
because we have used the recurrence relation to truncate the hierarchy, while Shoji and 
Komatsu have used ^3 = 0. The evolution equation for the shear can be further simplified 
by using Eq. (|3.3|) , leading to: 



a 



+ 



1 S 
3a 



G + ak 2 { 8 



l + w 



3c: 



(3.10) 



3.4 Estimating higher order momenta 



One way to close the system governing the fluid approximation is to replace G and £ by the 
usual quantities 6 and a multiplied by functions depending only on background quantities 



3 The recurrence relation in the massless limit is better motivated theoretically, since ^ oc ji (kr) when 
metric perturbations vanish or satisfy a simple constraint (namely, 4> + i[) = in the Newtonian gauge). In 
the massive case, the formal solution involves more complicated oscillating functions with arguments going 



from ~ kr in the massless limit to ~ (kr) 1 in the massive limit, as can be checked from eq. (2.£ 
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(in the same way that we already approximated 5p by c 2 5p). More explicitly, our aim is 
to write an approximation of the type S = 3w a a, where w a could be any function of time 
going from one third in the relativistic limit to zero in the non-relativistic one. Since 6 
and are not gauge- independent, we should search for a similar approximation holding on 
their gauge-independent counterpart. In the synchronous gauge, such an approximation 
would read 



@ + ak 2 [8— 3c?) = 3w e \6 + ak 2 ] . (3.11) 

\ 1 + tu y J J L J 



However, we will stick to the notations of ||], who introduced a viscosity speed related to 
our we through 

cl is = ^w e (l+w) . (3.12) 
With such assumptions, the approximate equation for the shear would read 



a = -3 ( - + - 

r a 



2 2 
~-c g -w 



a. 



a + xi^ \ 2 + 2ak 2 ] . (3.13) 

3 1 + w 



Since the suppression factor q 2 /e 2 which appears in Eq. ( |3.8| , 3.9) compared to Eq. ( |2.2c 



2.2d| ) is also found in the pressure integral compared to the energy density integral, we 



may guess that the relative behaviour is similar, i.e. related by w. This leads to a guess 
w a = w and wg = w which implies c 2 is = |u> (l + w). However, the same logic would imply 
c cff = w > which we have shown in Fig. |l] is not exactly true. 

Let us investigate a bit how to approximate higher momenta quantities like and E. 
If we want to approximate for instance, we may assume some functional form of ^i(q) 
described by a single (time dependent) parameter. We can make the ansatz ^f\(q/e) ~ 
a\ n (t) (f)™, and then use to determine the parameter a\ n (t). We then find 

6afl V n n S" ' (3 - 14) 

Jo /oW)d^(f) 

The guess c 2 is = jw (l + w) can be seen to be a special case of this approach having 
n = — 1. The problem is that the value of n best approximating the behavior of and 
other momenta is not the same in the relativistic and non-relativistic limit. In fact our 
testing shows that this guess sources a too much during the relativistic to non-relativistic 
transition compared to the exact solution. Instead we got much better results by using 
c vis = 3u>c^, which avoids this excessive sourcing during the transition, while still reducing 
to 1/3 in the relativistic limit. 

For the ratio E/er, the assumption of a ^-independent ^2 (i-e. n = 0) yields w a = 
p/(3p), which provides satisfactory results and is adopted in the schemes described below. 

We speculate that by pushing these kinds of considerations further, one could find 
better approximations for c^ s , c 2 - and w a . It is also possible that another independent 
equation could be found, and that it would allow a better determination of c e g. 
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3.5 Implementation in CLASS 

For comparison, we have implemented 3 different Non-Cold Dark Matter Fluid Approxi- 
mations (NCDMFA) in CLASS which differ only in their respective equation for the shear. 
In correspondence with the Ultra-relativistic Fluid Approximation discussed in H , we have 
named the approximations MB, Hu and CLASS: in the relativistic limit, they reduce to their 
relativistic counterpart in ||. In all three approximations we are using Eq. ( 3.2a| ) and ( p. 2b ) 
as the first two equations with c e ff = c„. The respective equations for the shear read 



0"Hu 



O" CLASS 



t a 



c„ 



ll 

3p 



cr + - 



4 ci 



31 + w; 



2(9 + /i + 6r) 



3wc n 



aw 3 1 + w 



20 + h + Qi] 



1 a 
- + - 

t a 



IP 

3p 



3 1 + w L 



29 + h 



w, 



3wc„ 



(3.15a) 

(3.15b) 
(3.15c) 



The second shear equation, named Hu, corresponds exactly to the prescription of Ref. || for 
approximating massive neutrinos. The first shear equation, MB, comes directly from ( 3. 13j) 
with the values of w a and c v j s motivated in the previous subsection. Finally, in 0, we found 
that removing the f] term leads to slightly better results for the matter power spectrum, 
and can be justified using an analytic approximation to the exact equations. By analogy, 
we also define in the massive neutrino case a CLASS approximation identical to the MB one 
except for the omission of this term. 

In Fig. [2] we have tested these three fluid approximations in a model with no massless 
neutrinos and 3 degenerate massive neutrinos. The three approximations work very well 
as long as the neutrinos are light and become non-relativistic after photon decoupling. 
Like in the massless case, the CLASS approximation is slightly better for predicting the 
matter power spectrum on small scales, and we set it to be the default method in the code. 
When the mass increases, the fluid approximation alters the CMB spectra on small angular 
scales (I > 2500), but the error remains tiny (only 0.02% for I = 2750 for three species with 
m = leV). The effect on the matter power spectrum is stronger: with three 1 eV neutrinos, 
the P(k) is wrong by 1 to 3% for k £ [0.05; l]fcMpc -1 . Hence, we recommend to use the fluid 
approximation for any value of the mass when computing CMB anisotropics, and only below 
a total mass of one or two eV's when computing the matter power spectrum. However, 
cosmological bounds on neutrino masses strongly disfavour larger values of the total mass. 
This means that in most projects, CLASS users can safely use the fluid approximation for 
fitting both CMB and large scale structure data. 



4. Standard massive neutrinos 

We first illustrate our approach with the simple case of standard massive neutrinos with 
a Fermi-Dirac distribution. In this case, for each neutrino, the user should provide two 
numbers in the input file: the mass m, and the relative temperature T_ncdm = T v /T^ 
(the ratio of neutrino to the photon temperature). The CLASS input file explanatory.ini 
recommends to use the value T_ncdm=0. 71599, which is "fudged" in order to provide a 
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Figure 2: On the left, we have shown the percentage difference in the Cj for three degenerate 
neutrino species with mass m = O.OOleV, m — O.OleV, m = O.leV and m = leV respectively, in 
runs with/ without the fluid approximation. The fluid approximation works very well as long as the 
neutrinos are relativistic, so this is what we expect. On the right we have shown the matter power 
spectrum for the same masses. Here the agreement is not so good as the mass becomes higher. 



mass-to-density ratio m/uj u = 93.14 eV in the non-relativistic limit. This number gives a 
very good approximation to the actual relic density of active neutrinos, resulting from an 
accurate study of neutrino decoupling [|llj. However, when comparing the CLASS results 
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with those from CAMB, we take T_ncdm=0.7133 in order to recover the mass-to-density ratio 
assumed in that code. Finally, if no temperature is entered, the code will default to the 
instantaneous decoupling value of (4/H) 1 ^ 3 . 

4.1 Agreement with CAMB 



o 




Figure 3: Relative difference between CAMB and CLASS spectra in a model with J7„ 
massless neutrinos, and reference accuracy settings. The two codes agree rather well. 



0.02, two 



In Fig. |3|, we compare the CMB and matter power spectrum from CAMB and CLASS 
(without the NCDM fluid approximation) for two massless and one massive neutrino with 
£l u = 0.02 (corresponding to a mass m — 0.923eV). We used high accuracy settings for 
CAMB, described in || under the name [CAMB:07]. For CLASS, we used the input file 
cl_ref .pre, which corresponds to the setting [CLASS:01J in || for parameters not related 
to NCDM; for the latter, cl_ref . pre contains the settings described in the first column 
of Table [l|. For such settings and in absence of massive neutrinos, the two temperature 
spectra would agree at the 0.01% level in the range I G [20; 3000]; at the 0.02% level for 
polarization in the same range; and at the 0.01% level for the matter power spectrum for 
k < 1/iMpc . With a neutrino mass close to 1 eV, we see in Fig. || that the discrepancy is 
approximately six times larger than in the massless case. However, it remains very small: 
even with massive neutrinos the two codes agree to better than 0.1% for the CMB and 
matter power spectra. This is by far sufficient for practical applications. 

In a perfect implementation of massless and massive neutrinos in Boltzmann codes, we 
expect that in the relativistic limit m <T„ (where T® is the neutrino temperature today) 
the spectra would tend towards those obtained with three massless species (provided that 
we are careful enough to keep the same number of relativistic degrees of freedom N e g). We 
performed this exercise for both codes, and the results are presented in Fig. |j. It appears 
that with a small enough mass, CLASS can get arbitrarily close to the fully relativistic case: 
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Figure 4: This is a test of how well CAMB and CLASS recovers the massless limit. We compute a 
model with S!„ = 1.2 • 10~ 4 and 3 massive neutrinos with degenerate mass. This setting corresponds 
to a neutrino mass of = 2.8 • 10 _4 eV, which is not exactly massless, but it is the best we 
can do since the mass parameter can not be set directly in CAMB. Setting the mass parameter in 
CLASS to rrii = 10~ 8 eV reveals that we are in part seeing the effect of the neutrino going slightly 
non-relativistic at late times. 



with a mass of 10 8 eV, the difference is at most of 0.03% in the Cj's and 0.05% in the 
P(k). This test is another way to validate the accuracy of our implementation. 

4.2 Accuracy settings 

We now come to the question of defining degraded accuracy settings for computing the 
spectra in a fast way, while keeping the accuracy of the results under control. For such 
an exercise, we need to define a measure a precision. Like in Q, we will use an effec- 
tive x 2 which mimics the sensitivity of a CMB experiment like Planck to temperature 
and E-polarisation anisotropies. Taking the runs with accuracy settings cl_ref .pre as 
a reference, we decrease the precision for each parameter while keeping the Ax 2 roughly 
below a given limit, chosen to be either 0.1 or 1. This exercise was already performed 
in H for all parameters not related to NCDM, leading to the definition of two precision 
files chi2pl0. 1 .pre and chi2pll.pre which are available on the CLASS web site. Here, 
we only need to set the NCDM precision parameters in these two files to correct values. 
Our results are listed in Table [l], in the second and third columns. They take advantage 
of the fluid approximation, and use an extremely small number of momenta (8 or 5 only). 
We checked that these settings provide the correct order of magnitude for Ax 2 within a 
wide range of neutrino masses, at least up to 2 eV. This is shown in Table [2] for the two 
cases chi2pl0 . 1 .pre and chi2pll .pre, as well as for the case chi2pll .pre with the fluid 
approximation removed. Around m = 2 eV, the error induced by the fluid approximation 
starts increasing significantly: when exploring this region, the user should either turn off 
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the approximation, or increase the value of the kr trigger. Given current limits on active 
neutrino masses, the interesting mass range to explore is below 2 eV, and in most projects, 
the CLASS users can safely employ the default settings of chi2pl0 . 1 .pre and chi2pll .pre 
including the fluid approximation. 

These settings are optimised for fitting the CMB spectra only. For the matter power 
spectra, the files chi2pl0. 1 .pre and chi2pll.pre produce an error of the order of a 
few per cents in the range k G [0.05; lJfoMpc" 1 (for any neutrino mass and with/without 
the fluid approximation). In order to get accurate matter power spectra, it is better 
to employ the settings cl_permille .pre, cl_2permille .pre, cl_3permille .pre, which 
lead to a precision of 1, 2 or 3 per mille for Cf T in the range 2 < I < 3000, even in 
the presence of neutrino masses. In these files, we fixed the fluid approximation trigger 
to a rather larger value in order to get a precision of one permille for the matter power 
spectrum for k < 0.2/iMpc -1 and m < 2 eV, or a bit worse for mildly non-linear scales 
k G [0.2;l]feMpc _1 . The power spectrum accuracy with such settings is indicated in Table H 
for various values of the mass. 





cl_ref.pre chi2pl0 . 1 .pre 


chi2pll .pre 


tol_ncdm_bg 


10 -io 


10" 5 


10" 5 


tol_ncdm 


10 -io 


10" 4 


10~ 3 


l_max_ncdm 


51 


16 


12 


fluid approximation 


none 


ncdmf a_class 


ncdmf a_class 


kr trigger 




30 


16 


number of q (back.) 


28 


11 


11 


number of q (pert.) 


28 


8 


5 


number of neutrino equations 


1428 


136^3 


65^3 



Table 1: Accuracy parameters related to NCDM in the three precision files cl_ref.pre, 
chi2pl0. 1 .pre and chi2pll .pre. When the fluid approximation is used, the method described in 
section |^ is employed, and the switching time is set by the above values of kr. Below these param- 
eters, we indicate the corresponding number of momenta sampled in background quantities and in 
perturbation quantities, as well as the number of neutrino perturbation equations integrated over 
time, equal to (l_max_ncdm+l) times the number of sampled momenta when the fluid approximation 
is not used, and to three afterwards. 

4.3 Performance 

The quadrature method reveals to be extremely useful since even with five values of the 
momenta, we get accurate results leading to 0.2%-0.3% accuracy on the C[s, 0.1% accuracy 
on the P{k) and Ax 2 ~ 1. Traditional Boltzmann codes employ 14 momenta in order to 
achieve a comparable precision. In the presence of massive neutrinos, the total execution 
time of a Boltzmann code is dominated by the integration of the perturbation equations, 
which depends on the total number of perturbed variables, itself dominated by the number 
of massive neutrino equations. By reducing the number of momenta from 14 to 5, the 
quadrature method speeds up the code by more than a factor two. We find that the use of 
the fluid approximation leads to an additional 25% speed up for standard accuracy settings 
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mass (pV) 


chi2pl0 . 1 . pre 


ch.i2pll . prG 


same without a/pprox. 


i(r 3 


0.087 


0.94 


0.90 


10~ 2 


0.087 


0.93 


0.92 


O.l 


0.092 


0.90 


0.92 


i 


0.083 


0.96 


0.82 


2 


0.157 


1.10 


0.93 



Table 2: For a CMB instrument with the sensitivity of Planck, \ 2 difference between the spectra 
obtained with reference accuracy settings and with degraded accuracy settings, for various values 
of the neutrino mass (all models have two massless and one massive neutrinos). This shows that 
our accuracy settings chi2pl0 . 1 .pre and chi2pll.pre always lead to an accuracy of roughly 
Ax 2 ~ 0.1 or Ax 2 ~ 1 respectively. The last column correspond to the settings of chi2pll .pre, 
but without the fluid approximation. 



mass (eV) 


k < 0.2/iMpc" 1 


k e [0.2;l]fcMpc -1 


10~ 3 


0.04% 


0.12% 


io- 2 


0.04% 


0.12% 


0.1 


0.05% 


0.12% 


1 


0.06% 


0.8% 


2 


0.2% 


1.5% 



Table 3: Maximum error induced by any of the cl_permille .pre, cl_2permille .pre or 
cl_3permille .pre precision settings on the linear matter power spectrum P(k), for approximately 
linear scales k < 0.2/i.Mpc (first column) or mildly non-linear scales k <G [0.2; l]ft,Mpc _1 (second 
column) , and for various values of the neutrino mass (all models have two massless and one massive 
neutrinos). The fluid approximation introduces an error which remains below the per mille level 
until k = O^/iMpc -1 for m < 2 eV, and exceeds this level for larger masses. 

(like those in the file chi2pll .pre). In total, for a single massive neutrino, our method 
speeds up the code by a factor 3. This means that instead of being 4.5 times slower in 
presence of one massive neutrino, CLASS only becomes 1.5 times slower. We checked these 
numbers with various masses and accuracy settings. 

4.4 Realistic mass schemes 

We have proved in this section that CLASS can be employed in any project requiring high- 
precision computations of cosmological observables in presence of massive neutrinos. It is 
of course perfectly suited for realistic situations with different neutrino species and masses. 

To illustrate this, we display in Figure |5| the ratio of pairs of matter power spectra for 
models with three massive neutrinos satisfying constraints from atmospheric/solar oscilla- 
tion experiments (Am^ = 7.6 x 10~ 5 eV 2 , Am| 2 = ±2.4 x 10~ 3 eV 2 ). Each pair of 
models corresponds to one normal hierarchy and one inverted hierarchy scenario, with the 
same total mass M u , equal to 0.100 eV, 0.115 eV or 0.130 eV. The first total mass is very 
close to the minimum allowed value for the inverted hierarchy, M u ~ 0.0994 eV. For each 
pair or models with a given M u : 
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Figure 5: Ratio of matter power spectra for pairs of models with three massive neutrinos, obeying 
either to the normal or inverted hierarchy scenario, but with a common total mass for each pair: 
M v = 0.100 eV, 0.115 eV or 0.130 eV. The various effects observed here are discussed in the text. 



on intermediate scales, the bump reflects the difference in the three free-streaming 
scales involved in the two models. 

in the large k limit, the two spectra are offset by 0.03% to 0.22%: it is known that in 
this limit, the suppression in the power spectrum induced by neutrino free-streaming 
depends mainly on the total mass (through the famous —8f u approximate formula), 



but also slightly on the mass splitting (in [ 13 1 , a more accurate formula gives the 
suppression as a function of both the total mass and number of degenerate massive 
neutrinos). When M u increases, the two models are less different from each other 
(they go towards a common limit, namely the degenerate mass scenario), and the 
discrepancy is less pronounced. 
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• in the small k limit, the two spectra are nearly identical. The tiny difference, which 
increases when M v decreases, is due to the fact that in the inverted hierarchy model, 
there is a very light neutrino just finishing to complete its non-relativistic transi- 
tion today. It therefore has a non-negligible pressure, which slightly affects metric 
perturbations on large wavelengths. 

Observing the difference between these two models would be extremely challenging, al- 



though 21 cm surveys could reach enough sensitivity [14]. 



5. Beyond standard massive neutrinos 

In this section we will illustrate the power and flexibility of the non-cold Dark Matter im- 
plementation in CLASS, by implementing different models which have already been studied 
elsewhere in the literature. 



5.1 Massive neutrinos with large non-thermal corrections 

It is plausible that some new physics can introduce non-thermal corrections to an otherwise 
thermal Fermi-Dirac distribution function. One might think of using CMB and large scale 



structure data to put bounds on such non-thermal corrections, as was described e.g. in [15]. 



CLASS is ideally suited for playing with such models. As a test case, we take the following 
distribution from 15]: 



Att 2 ( (q-q c 



\2 



(5.1) 



which is the Fermi-Dirac distribution with an added Gaussian peak in the number density. 
This distribution could presumably be the result of some particle suddenly decaying into 
neutrinos at a late time. 

In practise, we only need to change the expression for /(g) in CLASS, which appears in a 
unique line (in the function background_ncdm_distribution() ). All the rest, like density- 
to-mass relation and computation of the logarithmic derivative, is done automatically by 
the code. In particular, we do not need to change the accuracy parameters tol_ncdm and 
tol_ncdm_bg: the momentum sampling algorithm automatically increases the number of 
momenta by a significant amount, in order to keep the same precision. If this was not the 
case, the effect of the peak would be underestimated because of under sampling, and the 
parameter extraction would then likely be biased. 

In Fig. ^, we show the CMB and matter power spectra for this model, relative to a 
standard model with three thermally distributed neutrinos. The two models are chosen to 
share exactly the same masses and the same initial number of relativistic degrees of freedom 
N e g. Nevertheless, they do not have the same non-relativistic neutrino density and average 
neutrino momentum; in particular, non-thermal neutrinos in the decay peak become non- 
relativistic slightly later. This induces a combination of background and perturbation 
effects affecting CMB and matter power spectra in a significant way. 
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Figure 6: C;'s and P(k)'s for a model of 3 degenerate neutrinos with the non-thermal distribu- 
tion (5.1) using parameters m — 1.0 eV, A — 0.018, a = 1.0 and q c = 10.5. This corresponds to 



iVeff = 3.98486. We have compared this model to a model with degenerate thermal neutrinos with 
the same mass and N c ff. The signal is due to a combination of background and perturbation effects: 
although the mass and the relativistic density are the same, the non-relativistic density and the 
average momentum differ significantly in the two models. 



5.2 Warm dark matter with thermal-like distribution 

There is an infinity of possible warm dark matter models, since the phase-space distribution 
of warm dark matter depend on the details of its production mechanism. The most widely 
studied model is that of non-resonantly produced warm dark matter with a rescaled Fermi- 
Dirac distribution, having the same temperature as that of active neutrinos. This model 
is implemented in the default CLASS version: when the user enters a temperature, a mass 
and a density O nC( j m (or w nc d m ) for the same species, the code knows that the degeneracy 
parameter in front of the Fermi-Dirac distribution must be rescaled in order to match these 
three constraints simultaneously. The code will also ensure that the perturbations begin 
to be integrated when the non-cold species is still relativistic, in order to properly follow 
the transition to the non-relativistic regime. 

We illustrate this by running a AWDM model with a mass of m = IkeV or m = lOkeV 
and a density rincdm = 0.25, with or without the fluid approximation. We compare the 
results with those of ACDM with Slcdm = 0.25, in order to show the well-known suppression 
effect of WDM in the small-scale limit of the matter power spectrum. It appears that the 
fluid approximation works very well in those cases, unless one wants to resolve the details 



of the WDM acoustic oscillations on very small scales, first predicted in [16] 



5.3 Warm dark matter with non-trivial production mechanism 

Non-resonantly produced warm dark matter candidates are severely constrained by Lyman- 
a bounds, but such bounds do not apply to other warm particles which could have been 
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Figure 7: P(k)'s for a Warm Dark Matter model with m = 1 keV (left) and m = 10 keV (right). 
The fluid approximation can be seen to be a very good approximation in this case, though it does 
not catch the acoustic oscillations precisely. 



produced through more complicated mechanisms (e.g. resonant production), leading to a 
non-trivial, model-dependent phase-space distribution function fi7| . It is not always easy 
to find a good analytic approximation for such a distribution; this is anyway not an issue 
for CLASS, since the code can read tabulated values of f(p) from an input file. 

0.45 | 1 1 1 1 1 1 r-^ — i — h — i — 1 — h — i — — ~n 




q = p/T k (h/Mpc) 

Figure 8: P(k)'s for a Warm Dark Matter model with a non-trivial production mecha- 
nism for a mass of m = 2 keV compared to the same model with Cold Dark Matter. Note 
that the normalisation of the distribution function is arbitrary; when both m_ncdm and one of 
{0mega_ncdm, omega_ncdm} is present for some species, CLASS will normalise the distribution con- 
sistently. 
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We illustrate this case by taking a particular model for resonantly produced sterile 
neutrinos, which distribution was computed numerically by |l8| (simulating the details of 
sterile neutrino production and freeze-out), and stored in a file with discrete qi, fi values. 
Again, we only need to specify the name of this file in the CLASS input file, to enter a value 
for the mass and for the density ^ ncc i m , and the rest is done automatically by the code 
(finding the mass-density relation and the correct normalization factor for f(q), defining 
the new momentum steps, deriving [d In /] / [d In q] with a good enough accuracy). The 
assumed f{q) and the resulting matter power spectrum when the mass is set to m = 2keV 
is shown in Fig. [|. By eye, this spectrum seems identical to a thermal-like WDM one, but 
the cut-off is in fact much smoother due to an excess of low-momentum particles in this 
model (which behave like a small cold dark matter fraction). 

6. Conclusions 

A large fraction of the activity in cosmology consists in deriving bounds on particle physics 
in general, and on the neutrino and dark matter sector in particular. Fitting cosmological 
data with non-standard neutrinos or other non-cold relics require non-trivial changes in 
existing public Boltzmann codes. Moreover, running parameter extraction codes including 
massive neutrinos or more exotic non-cold relics is computationally expensive due to a 
significant increase in the number of differential equations to be solved numerically for 
each set of cosmological parameters. 

The newly released Cosmic Linear Anisotropy Solving System aims at rendering this 
task easy and fast. The code provides a very friendly and flexible input file in which users 
can specify a lot of non-standard properties for the NCDM sector: masses, temperatures, 
chemical potentials, degeneracy parameters, etc. Moreover, the Fermi-Dirac distribution 
function is not hard-coded in CLASS; it is just a default choice appearing in one line of 
the code, which can be very easily modified. Even when a non-thermal distribution f(q) 
does not have a simple analytic expression, the code can be told to read it directly from a 
file. After reading this function, CLASS performs a series of steps in a fully automatic way: 
finding the mass-density relation, defining an optimal sampling in momentum space with a 
sophisticated but fast algorithm, and accurately computing the derivative of f(q), needed 
in the perturbation equations. 

In this paper, we presented the main two improvements related to the NCDM sector 
in CLASS: an adaptive quadrature sampling algorithm, which is useful both for the purpose 
of flexibility (the sampling is always adapted to any new distribution function) and speed 
(the code sticks to a minimum number of momenta, and hence, of perturbation equations); 
and a fluid approximation switched on inside the Hubble radius. We showed that the latter 
approximation works very well for realistic active neutrinos (with a total mass smaller than 
1 — 2eV), and for warm dark matter candidates becoming non-relativistic during radiation 
domination. In between these two limits, there is a range in which the accuracy of the 
fluid approximation is not well established, and in which the user may need to keep the 
approximation off, at the expense of increasing the execution time. However, the range 
between a few eV and few keV is usually not relevant in most realistic scenarios. 
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The adaptive quadrature sampling algorithm and the fluid approximation both con- 
tribute to a reduction in the total execution time of the code by a factor of three for 
ordinary neutrinos. This means that when one massive neutrino species is added to the 
ACDM model, CLASS becomes 1.5 times slower instead of 4.5 times slower like other codes. 
Since the code is already quite fast in the the massless case, we conclude that the global 
speed up is significant and appreciable when fitting cosmological data. 
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